# ============================================================================
# 非线性GPU性能测试 - 材料非线性+几何非线性
# 测试GPU矩阵组装在非线性分析中的性能，包含刚度矩阵实时更新
# ============================================================================

wipe
puts "非线性GPU vs CPU性能对比测试 - 刚度矩阵实时更新"

# ============================================================================
# 非线性测试模型配置
# ============================================================================

model basic -ndm 2 -ndf 2

# 优化网格设计：减少条件数问题
set nEleX 50
set nEleY 50
set dx 0.1
set dy 0.1

puts "创建非线性测试模型: ${nEleX}x${nEleY} = [expr $nEleX*$nEleY] 个单元"

# 创建节点
set nodeTag 1
for {set j 0} {$j <= $nEleY} {incr j} {
    for {set i 0} {$i <= $nEleX} {incr i} {
        node $nodeTag [expr $i * $dx] [expr $j * $dy]
        incr nodeTag
    }
}

# ============================================================================
# 选择监测节点用于CPU vs GPU结果对比
# ============================================================================

# 选择顶部中心节点作为监测点
set monitorNode [expr $nEleY*($nEleX+1) + ($nEleX+1)/2]
puts "监测节点: $monitorNode (顶部中心节点，坐标: [expr ($nEleX+1)/2 * $dx], [expr $nEleY * $dy])"
puts "此节点的响应将用于CPU与GPU计算结果精度对比"

# ============================================================================
# 非线性材料模型 - Drucker-Prager + 大变形
# ============================================================================

# 非线性土体材料 - DruckerPrager模型
# 参数: K, G, sigma_y, rho, rhoBar, Kinf, Ko, delta1, delta2, H, theta
nDMaterial DruckerPrager 1 2 \
    100.0e6  75.0e6  50.0e3 \
    0.05 0.1 0.0 0.0 \
    0.0 0.0 0.0 1.0 \
    1800.0

puts "材料: DruckerPrager非线性土体模型 (会产生材料非线性)"

# 创建非线性四边形单元
set eleTag 1
for {set j 0} {$j < $nEleY} {incr j} {
    for {set i 0} {$i < $nEleX} {incr i} {
        set n1 [expr $j*($nEleX+1) + $i + 1]
        set n2 [expr $j*($nEleX+1) + $i + 2]
        set n3 [expr ($j+1)*($nEleX+1) + $i + 2]
        set n4 [expr ($j+1)*($nEleX+1) + $i + 1]
        
        # 使用quad单元，支持材料非线性
        element quad $eleTag $n1 $n2 $n3 $n4 1.0 "PlaneStrain" 1
        incr eleTag
    }
}

# ============================================================================
# 边界条件和加载
# ============================================================================

# 底部固定边界
for {set i 1} {$i <= [expr $nEleX+1]} {incr i} {
    fix $i 1 1
}

# 侧边边界：允许竖向位移，限制水平位移  
for {set j 1} {$j <= $nEleY} {incr j} {
    set leftNode [expr $j*($nEleX+1) + 1]
    set rightNode [expr $j*($nEleX+1) + ($nEleX+1)]
    fix $leftNode 1 0
    fix $rightNode 1 0
}

# 非线性加载模式 - 大载荷产生材料非线性
pattern Plain 1 Linear {
    # 顶部施加大的分布荷载，确保进入塑性
    for {set i 1} {$i <= [expr $nEleX+1]} {incr i} {
        set topNode [expr $nEleY*($nEleX+1) + $i]
        # 大荷载：1000 kN/m²，确保材料进入非线性状态
        load $topNode 0.0 -10000.0
    }
}

set totalDOF [expr ($nEleX+1)*($nEleY+1)*2]
puts "模型: [expr $nEleX*$nEleY] 单元, ${totalDOF} 自由度"
puts "预期：材料非线性导致刚度矩阵每步更新"

# ============================================================================
# CPU 非线性分析
# ============================================================================

puts "\n=== CPU 非线性性能测试 ==="

constraints Plain
numberer RCM  ; # 使用RCM重新编号减少带宽
system BandGeneral
test NormDispIncr 1.0e-6 50  ; # 增加最大迭代次数
algorithm Newton
integrator LoadControl 0.05  ; # 更小的载荷步
analysis Static

# ============================================================================
# CPU结果记录器设置
# ============================================================================

# 删除已存在的输出文件
file delete -force "cpu_response.txt"

# 创建记录器记录监测节点的位移响应
recorder Node -file "cpu_response.txt" -time -node $monitorNode -dof 1 2 disp

puts "CPU分析: 监测节点 $monitorNode 的响应将记录到 cpu_response.txt"
puts "记录内容: 时间, X位移, Y位移"

set cpuTotalTime 0
set testSteps 20  ; # 更多步数测试非线性

puts "开始CPU非线性分析 ($testSteps 步)..."
puts "注意：每步都会重新组装刚度矩阵"

for {set i 1} {$i <= $testSteps} {incr i} {
    set stepStart [clock clicks -milliseconds]
    set ok [analyze 1]
    set stepEnd [clock clicks -milliseconds]
    set stepTime [expr $stepEnd - $stepStart]
    set cpuTotalTime [expr $cpuTotalTime + $stepTime]
    
    if {$i % 5 == 0} {
        puts "  CPU step $i: ${stepTime} ms (非线性分析)"
    }
    
    if {$ok != 0} {
        puts "CPU非线性分析收敛困难，在第 $i 步"
        # 尝试减小载荷步长
        integrator LoadControl 0.01
        set ok [analyze 1]
        if {$ok != 0} {
            puts "CPU分析失败在第 $i 步"
            break
        }
        # 恢复载荷步长
        integrator LoadControl 0.05
    }
}

set cpuAvgTime [expr $cpuTotalTime / $testSteps]
puts "CPU非线性结果: 总时间=${cpuTotalTime}ms, 平均=${cpuAvgTime}ms/步"

# 完全清理并重建模型以确保GPU分析从初始状态开始
wipe

# ============================================================================
# 重建模型用于GPU分析
# ============================================================================

puts "\n重建模型用于GPU分析..."
model basic -ndm 2 -ndf 2

# 重新创建节点
set nodeTag 1
for {set j 0} {$j <= $nEleY} {incr j} {
    for {set i 0} {$i <= $nEleX} {incr i} {
        node $nodeTag [expr $i * $dx] [expr $j * $dy]
        incr nodeTag
    }
}

# 重新定义材料
nDMaterial DruckerPrager 1 2 \
    100.0e6  75.0e6  50.0e3 \
    0.05 0.1 0.0 0.0 \
    0.0 0.0 0.0 1.0 \
    1800.0

# 重新创建单元
set eleTag 1
for {set j 0} {$j < $nEleY} {incr j} {
    for {set i 0} {$i < $nEleX} {incr i} {
        set n1 [expr $j*($nEleX+1) + $i + 1]
        set n2 [expr $j*($nEleX+1) + $i + 2]
        set n3 [expr ($j+1)*($nEleX+1) + $i + 2]
        set n4 [expr ($j+1)*($nEleX+1) + $i + 1]
        element quad $eleTag $n1 $n2 $n3 $n4 1.0 "PlaneStrain" 1
        incr eleTag
    }
}

# 重新施加边界条件
for {set i 1} {$i <= [expr $nEleX+1]} {incr i} {
    fix $i 1 1
}
for {set j 1} {$j <= $nEleY} {incr j} {
    set leftNode [expr $j*($nEleX+1) + 1]
    set rightNode [expr $j*($nEleX+1) + ($nEleX+1)]
    fix $leftNode 1 0
    fix $rightNode 1 0
}

# 重新定义荷载模式
pattern Plain 1 Linear {
    for {set i 1} {$i <= [expr $nEleX+1]} {incr i} {
        set topNode [expr $nEleY*($nEleX+1) + $i]
        load $topNode 0.0 -10000.0
    }
}

puts "模型重建完成，GPU分析将从初始状态开始"

# ============================================================================
# GPU 非线性分析 - 测试GPU矩阵组装更新性能
# ============================================================================

puts "\n=== GPU 非线性性能测试 ==="
puts "GPU矩阵组装将在每个Newton迭代中更新刚度矩阵"

constraints Plain  
numberer RCM
system GPUSolver  ; # 使用GPU求解器
test NormDispIncr 1.0e-6 50
algorithm Newton  
integrator LoadControl 0.05
analysis Static

# 重置时间为0以确保从初始状态开始
setTime 0.0

# ============================================================================
# GPU结果记录器设置
# ============================================================================

# 删除已存在的输出文件
file delete -force "gpu_response.txt"

# 创建记录器记录监测节点的位移响应
recorder Node -file "gpu_response.txt" -time -node $monitorNode -dof 1 2 disp

puts "GPU分析: 监测节点 $monitorNode 的响应将记录到 gpu_response.txt"
puts "记录内容: 时间, X位移, Y位移"

set gpuTotalTime 0

puts "开始GPU非线性分析 ($testSteps 步)..."
puts "监控GPU矩阵组装性能:"
puts ""

for {set i 1} {$i <= $testSteps} {incr i} {
    puts "--- GPU 非线性步 $i ---"
    set stepStart [clock clicks -milliseconds]
    set ok [analyze 1]
    set stepEnd [clock clicks -milliseconds]
    set stepTime [expr $stepEnd - $stepStart]
    set gpuTotalTime [expr $gpuTotalTime + $stepTime]
    
    puts "GPU步 $i: ${stepTime} ms (wall time, 包含刚度矩阵更新)"
    puts "^^ 上方应显示GPU矩阵组装详细耗时 ^^"
    puts ""
    
    if {$ok != 0} {
        puts "GPU非线性分析收敛困难，在第 $i 步"
        # 尝试减小载荷步长  
        integrator LoadControl 0.01
        set ok [analyze 1]
        if {$ok != 0} {
            puts "GPU分析失败在第 $i 步"
            break
        }
        # 恢复载荷步长
        integrator LoadControl 0.05
    }
}

set gpuAvgTime [expr $gpuTotalTime / $testSteps]
puts "GPU非线性结果: 总时间=${gpuTotalTime}ms, 平均=${gpuAvgTime}ms/步"

# ============================================================================
# 非线性性能分析和对比
# ============================================================================

puts "\n=========================================="  
puts "非线性GPU vs CPU性能对比结果"
puts "=========================================="
puts "测试模型: [expr $nEleX*$nEleY] 单元, ${totalDOF} 自由度"
puts "分析类型: 非线性静态分析 (材料非线性 + 刚度矩阵更新)"
puts "测试步数: $testSteps"
puts ""

puts "CPU 非线性性能:"
puts "  总时间: ${cpuTotalTime} ms" 
puts "  平均/步: ${cpuAvgTime} ms"
puts "  特点: 每步重新组装CPU刚度矩阵"
puts ""

puts "GPU 非线性性能:"
puts "  总时间: ${gpuTotalTime} ms"
puts "  平均/步: ${gpuAvgTime} ms"  
puts "  特点: 每步GPU并行组装刚度矩阵"
puts ""

if {$cpuTotalTime > 0 && $gpuTotalTime > 0} {
    set speedup [expr double($cpuTotalTime) / $gpuTotalTime]
    puts "非线性分析性能对比:"
    if {$speedup > 1.0} {
        puts "  ✓ GPU加速: [format %.2fx $speedup] 倍于CPU"
        puts "  ✓ 节省时间: [expr $cpuTotalTime - $gpuTotalTime] ms"
        puts "  ✓ GPU矩阵组装在非线性分析中显示优势"
    } else {
        set slowdown [expr double($gpuTotalTime) / $cpuTotalTime]  
        puts "  ⚠ GPU较慢: [format %.2fx $slowdown] 倍于CPU"
        puts "  ⚠ 额外时间: [expr $gpuTotalTime - $cpuTotalTime] ms"
    }
}

puts ""
puts "=========================================="
puts "非线性分析GPU矩阵组装分析"
puts "=========================================="
puts "在上方GPU测试输出中查看："
puts ""
puts "关键指标："
puts "• Matrix Assembly时间和百分比"
puts "• GPU Computation vs Memory Transfer比例"
puts "• 每步的Newton迭代次数" 
puts "• 每次Newton迭代的矩阵组装耗时"
puts ""
puts "非线性优势："
puts "• 刚度矩阵每步更新 → GPU并行组装显示优势"
puts "• Newton迭代多次求解 → CG求解器重用分解"
puts "• 大规模问题 → GPU内存带宽优势明显"
puts ""
puts "预期结果："
puts "• GPU矩阵组装应比线性分析显示更大优势"
puts "• 每步包含多次Newton迭代和矩阵更新"
puts "• 总体GPU性能应优于纯线性分析"

puts ""
puts "=========================================="
puts "非线性GPU性能测试完成!"
puts "检查上方GPU矩阵组装详细耗时分析"
puts "=========================================="

# ============================================================================
# CPU vs GPU 计算结果精度对比
# ============================================================================

puts ""
puts "=========================================="
puts "CPU vs GPU 计算结果精度对比分析"
puts "=========================================="

# 检查输出文件是否存在
if {[file exists "cpu_response.txt"] && [file exists "gpu_response.txt"]} {
    puts "比较监测节点 $monitorNode 的位移响应..."
    puts ""
    
    # 读取CPU结果
    set cpu_file [open "cpu_response.txt" r]
    set cpu_data [read $cpu_file]
    close $cpu_file
    
    # 读取GPU结果  
    set gpu_file [open "gpu_response.txt" r]
    set gpu_data [read $gpu_file]
    close $gpu_file
    
    # 分析数据行数
    set cpu_lines [split [string trim $cpu_data] "\n"]
    set gpu_lines [split [string trim $gpu_data] "\n"]
    
    puts "数据点数量:"
    puts "  CPU结果: [llength $cpu_lines] 个时间步"
    puts "  GPU结果: [llength $gpu_lines] 个时间步"
    
    if {[llength $cpu_lines] == [llength $gpu_lines] && [llength $cpu_lines] > 0} {
        puts ""
        puts "逐步位移对比 (前10步和最后5步):"
        puts "格式: 时间 | CPU_X位移 | GPU_X位移 | CPU_Y位移 | GPU_Y位移 | X差值 | Y差值"
        puts "-------------------------------------------------------------------"
        
        set max_diff_x 0.0
        set max_diff_y 0.0
        set total_diff_x 0.0
        set total_diff_y 0.0
        set compare_count 0
        
        # 对比前10步
        set show_steps [expr min(10, [llength $cpu_lines])]
        for {set i 0} {$i < $show_steps} {incr i} {
            set cpu_line [split [lindex $cpu_lines $i]]
            set gpu_line [split [lindex $gpu_lines $i]]
            
            if {[llength $cpu_line] >= 3 && [llength $gpu_line] >= 3} {
                set time_cpu [lindex $cpu_line 0]
                set time_gpu [lindex $gpu_line 0]
                set x_cpu [lindex $cpu_line 1]
                set x_gpu [lindex $gpu_line 1]
                set y_cpu [lindex $cpu_line 2]
                set y_gpu [lindex $gpu_line 2]
                
                set diff_x [expr abs($x_cpu - $x_gpu)]
                set diff_y [expr abs($y_cpu - $y_gpu)]
                
                puts [format "步%2d: %.3f | %12.6e | %12.6e | %12.6e | %12.6e | %.2e | %.2e" \
                      [expr $i+1] $time_cpu $x_cpu $x_gpu $y_cpu $y_gpu $diff_x $diff_y]
                
                set max_diff_x [expr max($max_diff_x, $diff_x)]
                set max_diff_y [expr max($max_diff_y, $diff_y)]
                set total_diff_x [expr $total_diff_x + $diff_x]
                set total_diff_y [expr $total_diff_y + $diff_y]
                incr compare_count
            }
        }
        
        # 如果总步数大于15，显示最后5步
        if {[llength $cpu_lines] > 15} {
            puts "  ..."
            set start_idx [expr [llength $cpu_lines] - 5]
            for {set i $start_idx} {$i < [llength $cpu_lines]} {incr i} {
                set cpu_line [split [lindex $cpu_lines $i]]
                set gpu_line [split [lindex $gpu_lines $i]]
                
                if {[llength $cpu_line] >= 3 && [llength $gpu_line] >= 3} {
                    set time_cpu [lindex $cpu_line 0]
                    set x_cpu [lindex $cpu_line 1]
                    set x_gpu [lindex $gpu_line 1]
                    set y_cpu [lindex $cpu_line 2]
                    set y_gpu [lindex $gpu_line 2]
                    
                    set diff_x [expr abs($x_cpu - $x_gpu)]
                    set diff_y [expr abs($y_cpu - $y_gpu)]
                    
                    puts [format "步%2d: %.3f | %12.6e | %12.6e | %12.6e | %12.6e | %.2e | %.2e" \
                          [expr $i+1] $time_cpu $x_cpu $x_gpu $y_cpu $y_gpu $diff_x $diff_y]
                    
                    set max_diff_x [expr max($max_diff_x, $diff_x)]
                    set max_diff_y [expr max($max_diff_y, $diff_y)]
                    set total_diff_x [expr $total_diff_x + $diff_x]
                    set total_diff_y [expr $total_diff_y + $diff_y]
                    incr compare_count
                }
            }
        }
        
        puts ""
        puts "精度对比统计:"
        if {$compare_count > 0} {
            set avg_diff_x [expr $total_diff_x / $compare_count]
            set avg_diff_y [expr $total_diff_y / $compare_count]
            puts [format "  X方向位移差值: 最大=%.6e, 平均=%.6e" $max_diff_x $avg_diff_x]
            puts [format "  Y方向位移差值: 最大=%.6e, 平均=%.6e" $max_diff_y $avg_diff_y]
            
            # 判断精度等级
            set max_diff [expr max($max_diff_x, $max_diff_y)]
            if {$max_diff < 1e-12} {
                puts "  ✓ 精度评级: 优秀 (差值 < 1e-12)"
            } elseif {$max_diff < 1e-10} {
                puts "  ✓ 精度评级: 很好 (差值 < 1e-10)"
            } elseif {$max_diff < 1e-8} {
                puts "  ✓ 精度评级: 良好 (差值 < 1e-8)"
            } elseif {$max_diff < 1e-6} {
                puts "  ⚠ 精度评级: 一般 (差值 < 1e-6)"
            } else {
                puts "  ❌ 精度评级: 较差 (差值 >= 1e-6)"
            }
            
            puts ""
            puts "结果文件已保存:"
            puts "  CPU结果: cpu_response.txt"
            puts "  GPU结果: gpu_response.txt"
            puts "  可用外部工具进一步分析差异"
        } else {
            puts "❌ 无法解析响应数据进行对比"
        }
    } else {
        puts "❌ CPU和GPU结果文件数据行数不匹配或文件为空"
        puts "  CPU: [llength $cpu_lines] 行"
        puts "  GPU: [llength $gpu_lines] 行"
    }
} else {
    puts "❌ 找不到结果文件进行对比"
    puts "  CPU文件存在: [file exists cpu_response.txt]"
    puts "  GPU文件存在: [file exists gpu_response.txt]"
    puts ""
    puts "请检查分析是否正常完成"
}

puts ""
puts "=========================================="
puts "测试总结"
puts "=========================================="
puts "1. 性能对比: 查看上方GPU vs CPU耗时结果"
puts "2. 精度验证: 查看上方CPU vs GPU数值对比"
puts "3. 结果文件: cpu_response.txt 和 gpu_response.txt"
puts "4. 监测节点: $monitorNode (顶部中心，最大变形位置)"
puts ""
puts "建议进一步分析:"
puts "• 使用绘图软件比较两个响应曲线"
puts "• 检查GPU调试输出中的向量操作性能"
puts "• 验证非线性收敛行为的一致性"